Geometry Optimization

This tutorial will demonstrate the use of PyBigDFT for performing geometry optimization. We will start with a simple example of a molecular system. Then we will move to an advanced example involving a slab of NaCl.

Molecule Example

Let’s begin with a simple example of an Aflatoxin B1 molecule.

cano = "CN1C=NC2=C1C(=O)N(C(=O)N2C)C"

We can use openbabel to create a system from its cannonical smiles representation.

from openbabel.openbabel import OBMol, OBConversion, OBBuilder
conv = OBConversion()
conv.SetInFormat("can")

mol = OBMol()
conv.ReadString(mol, cano)
mol.AddHydrogens()

builder = OBBuilder()
builder.Build(mol)

True

Then we convert from babel to PyBigDFT. An openbabel structure after calling build is usually not a good starting point. Normally you would do some combination of optimization and conformer search. For our case, we will do geometry optimization followed by a short molecular dynamics run. This will give us a decent geometry, but one that is not the actual minimum.

from BigDFT.Interop.BabelInterop import convert_babel_to_system, molecular_dynamics, optimize_system
sys_start = convert_babel_to_system(mol)
sys_opt = optimize_system(sys_start)
sys = molecular_dynamics(sys_opt, 10000, 300)

/Users/wddawson/Documents/CEA/binaries/bds/install/lib/python3.7/site-packages/BigDFT/IO.py:371: UserWarning: Unsupported bond type had to be set to 1 (i.e. aromatic)
  UserWarning)
/Users/wddawson/Documents/CEA/binaries/bds/install/lib/python3.7/site-packages/BigDFT/IO.py:371: UserWarning: Unsupported bond type had to be set to 1 (i.e. aromatic)
  UserWarning)

from BigDFT.Visualization import InlineVisualizer
viz = InlineVisualizer(400, 300)
viz.display_system(sys_start, colordict={x: "black" for x in sys_start}, show=False)
viz.display_system(sys, colordict={x: "blue" for x in sys}, show=True)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Now we are ready to optimize the geometry of this system suing BigDFT. First, set up a basic calculator and input file.

from BigDFT.Calculators import SystemCalculator
code = SystemCalculator(verbose=False, skip=True)

from BigDFT.Inputfiles import Inputfile
inp = Inputfile()
inp.set_xc("PBE")
inp.set_hgrid("0.37")

BigDFT has a number of built-in geometry optimization methods which we can probe from the documentation.

from BigDFT.InputActions import optimize_geometry
help(optimize_geometry)

Help on function optimize_geometry in module BigDFT.InputActions:

optimize_geometry(inp, method='FIRE', nsteps=50, betax=4.0, frac_fluct=1.0, forcemax=0)
    Optimize the geometry of the system

    Args:
       nsteps (int): maximum number of atomic steps.
       method (str): Geometry optimizer. Available keys:
          * SDCG:   A combination of Steepest Descent and Conjugate Gradient
          * VSSD:   Variable Stepsize Steepest Descent method
          * LBFGS:  Limited-memory BFGS
          * BFGS:   Broyden-Fletcher-Goldfarb-Shanno
          * PBFGS:  Same as BFGS with an initial Hessian obtained from a force
                    field
          * DIIS:   Direct inversion of iterative subspace
          * FIRE:   Fast Inertial Relaxation Engine as described by Bitzek et
                    al.
          * SBFGS:  SQNM minimizer, keyword deprecated, will be replaced by
                    SQNM in future release
          * SQNM:   Stabilized quasi-Newton minimzer
       betax (float): the step size for the optimization method.
          This stepsize is system dependent and it has therefore to be
          determined for each system.
       frac_fluct (float): Fraction of force fluctuations. Stop if
          fmax < forces_fluct*frac_fluct.
       forcemax (float): Max forces criterion when stop.

Let’s activate the SQNM method. The important parameter to adjust is betax which is the step length, and the appropriate value can be sensitive to the method and system being computed. BigDFT automatically computes a measure of the force fluctations that come from the discretization and SCF error. This is used as a convergence measure. If you want to converge tighter, you should consider decreasing hgrid and the wavefunction convergence (set_wavefunction_convergence)

inp.optimize_geometry(method="SQNM", betax=1.0)

And run.

log = code.run(input=inp, posinp=sys.get_posinp(), name="caf", run_dir="work")

Found 27 different runs

We will also do a calculation using implicit solvent.

from copy import deepcopy
inp2 = deepcopy(inp)
inp2.set_implicit_solvent()

log_imp = code.run(input=inp2, posinp=sys.get_posinp(), name="caf-imp", run_dir="work")

Found 55 different runs

The first thing we want to do is extract the energies at each step of the geometry optimization and plot them.

energy = []
forces = []
for run in log:
    energy.append(run.energy)
    forces.append(run.forcemax)

from matplotlib import pyplot as plt
fig, axs = plt.subplots(1,1)
axs2 = axs.twinx()
axs.plot(energy, 'kx', label="energy")
axs2.plot(forces, 'r+', label="force")
axs.set_xlabel("Iteration", fontsize=12)
axs.set_ylabel("Energy (Hartree)", fontsize=12)
axs2.set_ylabel("Force (A.U.)", fontsize=12)
axs.legend(loc="upper center")
axs2.legend(loc="upper right")
axs.ticklabel_format(useOffset=False)

We have to take some care with plotting the values from the calculation using implicit solvent. This is because sometimes BigDFT restarts a calculation if it senses that the extrapolated geometry guess is not working good. In that case, one of our logfile instances won’t have an energy attribute. A try catch can handle this smoothly.

energy = []
forces = []
for run in log_imp:
    try:
        energy.append(run.energy)
        forces.append(run.forcemax)
    except AttributeError as e:
        pass

from matplotlib import pyplot as plt
fig, axs = plt.subplots(1,1)
axs2 = axs.twinx()
axs.plot(energy, 'kx', label="energy")
axs2.plot(forces, 'r+', label="force")
axs.set_xlabel("Iteration", fontsize=12)
axs.set_ylabel("Energy (Hartree)", fontsize=12)
axs2.set_ylabel("Force (A.U.)", fontsize=12)
axs.legend(loc="upper center")
axs2.legend(loc="upper right")
axs.ticklabel_format(useOffset=False)

Next we will extract the geometry.

systems = []
for step in log:
    systems.append(deepcopy(sys))
    systems[-1].update_positions_from_dict(step.astruct)

systems_imp = []
for step in log_imp:
    systems_imp.append(deepcopy(sys))
    systems_imp[-1].update_positions_from_dict(step.astruct)

We can create a picture now with all three computed geometries overlapping.

viz = InlineVisualizer(400, 300)
viz.display_system(systems[0], colordict={x: "black" for x in systems[-1]}, show=False)
viz.display_system(systems[-1], colordict={x: "red" for x in systems[-1]}, show=False)
viz.display_system(systems_imp[-1], colordict={x: "blue" for x in systems[-1]}, show=True)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Slab Example

In this section, we will study a more sophisticated example using a slab of NaCl as our test case. Here will show how to do a constrained optimization by fixing the positions or our slab, and allowing a molecule to move freely on top of it. First, will create a slab of NaCl built with the helper of the Atomic Simulation Environment.

from ase.build import bulk

atoms = bulk('NaCl', 'rocksalt', a=5.64, orthorhombic=True)
atoms *= [4, 4, 4]

Now we convert to the BigDFT system format. Note that we set the middle cell value to infinity so that we can explore this system as a surface (instead of the bulk).

from BigDFT.Interop.ASEInterop import ase_to_bigdft
from BigDFT.Systems import System
from BigDFT.UnitCells import UnitCell
sys = System()
sys["SUR:1"] = ase_to_bigdft(atoms)
sys.cell = UnitCell([float(atoms.cell[0, 0]), float("inf"), float(atoms.cell[2, 2])], units="angstroem")

We won’t be optimizing the lattice constant, instead we want to optimize something that is sticking to the surface.

from BigDFT.IO import XYZReader
from BigDFT.Fragments import Fragment
with XYZReader("O2") as ifile:
    sys["ABS:2"] = Fragment(xyzfile=ifile)
sys["ABS:2"].translate([x - y for x, y in zip(sys["SUR:1"].centroid, sys["ABS:2"].centroid)])
sys["ABS:2"].translate([0, 0.4*sys.cell[2, 2], 0])

viz = InlineVisualizer(400, 300)
viz.display_system(sys)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol